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Extreme mass ratio binary systems, binaries involving stellar mass objects orbiting massive black 
holes, are considered to be a primary source of gravitational radiation to be detected by the space- 
based interferometer LISA. The numerical modelling of these binary systems is extremely challenging 
because the scales involved expand over several orders of magnitude. One needs to handle large 
wavelength scales comparable to the size of the massive black hole and, at the same time, to resolve 
the scales in the vicinity of the small companion where radiation reaction effects play a crucial role. 
Adaptive finite element methods, in which quantitative control of errors is achieved automatically 
by finite element mesh adaptivity based on posteriori error estimation, are a natural choice that 
has great potential for achieving the high level of adaptivity required in these simulations. To 
demonstrate this, we present the results of simulations of a toy model, consisting of a point-like 
source orbiting a black hole under the action of a scalar gravitational field. 
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I. INTRODUCTION 

As we enter the era of Gravitational Wave Astronomy, a number of experiments to detect and study gravitational 
waves have been set up and some others are presently under development (see Q for a recent account). Among the 
second group we find one of the experiments that ispresently attracting a considerable amount of attention: the 
Laser Interferometer Space Antenna (LISA) Q, |M ^ |^ , a collaboration between ES A and NASA that is scheduled 
to be launched in the next decade. Extreme-Mass-Ratio Binaries (EMRBs) are considered to be a primary source of 
gravitational radiation to be detected by LISA 6, 7]. They consist of a "small" object, such a main sequence star, 
a stellar mass black hole, or a neutron star, with mass m ranging from IMq to IO^Mq, orbiting a massive black 
hole (MBH) with mass M ranging from 10'^ (if we consider the case of intermediate mass black holes) to 10® Mq 
(the case of big supermassive black holes sitting in the center of galaxies). This translates to EMRBs with mass 
ratios, fi = m/M , in the range 10"'^ — 10~® . In order to exploit this type of systems through LISA, it is crucial to 
have a good theoretical understanding of the evolution of these systems, good enough to produce accurate waveform 
templates in support of data analysis efforts. 

Because there is no significant coupling between the strong curvature effects between the MBH and its companion, 
relativistic perturbation theory is a well suited tool to study EMBRs. Clearly, the accuracy of this approximation 
depends on the smallness of the ratio /i. The goal is to study the perturbations generated by the small body in 
the (background) gravitational field of the MBH, and how these perturbations affect the motion of the small body 
itself, which follows the geodesies of the perturbed spacetime. That is, one is after studying how the presence of 
the small body affects its own trajectory. This problem is usually known in literature as the radiation reaction 
problem. This is an old problem and several a ppr oaches to deal with it have been proposed (see the recent reviews 
by Poisson jlll , Detweiler IToll and Mino [I3,lil[l3). A pragmatic approach is to use energy-momentum balance 
arguments |l4l Il5l Il6l llTl llSl Il9||. Under this approach, one estimates the changes in the small body constants 
of motion by computing the fluxes of energy and angular momentum at infinity and through the MBH horizon. 
This approach works well in the adiabatic regime, when the time scale of the radiation reaction is much bigger 
than the orbital time scale. Until now it has dealt with special orbits of the Kerr black hole, but it has not yet 
produced results for generic orbits because of the difficulty of adjusting the third constant of motion, the Carter 
constant, present in generic geodesies of the Kerr spacetime. However, there have been some recent advances in this 
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direction !23, IH [H HI] . 

An alternative approach consists in trying to describe the radiation reaction effects on the small body as the 
action of a local self-force that is responsible for the deviations from the geodesic motion. A rigorous formulation 
of this concept has been given by the first time by Mino, Sasaki and Tanaka '24], and later, adopting an axiomatic 
approach, by Quinn and Wald (25.] . These works give a formal prescription to compute the self- force. For the practical 
implementation of these prescription some techniques have been proposed (see \2(^ for a recent progress report): the 
mode-sum scheme [2^ l28l l29l |3C| , and a regularization scheme based on zeta- function regularization techniques l3ll| . 

In this paper we want to advocate an alternative, and at the same time complementary, approach: The direct 
numerical integration of the (linearized) evolution equations for the gravitational field together with the evolution 
equations for the small body, which form a system of partial differential equations (PDEs) coupled to a system of 
ordinary differential equations (ODEs). We can already find in the literature an attempt to use numerical methods for 
simulating EMRBs |32| . In this work, the perturbations are computed by using the full general relativistic equations in 
the framework of the characteristic formulation (in contrast to a Cauchy-type formulation) , with the small body being 
described by an energy-momentum distribution that moves "rigidly" along a geodesic of the numerically computed 
spacetime. The main drawback of this work is the size of the small body, which is bigger than the MBH horizon. 
Another approach that has been used in the literature 33] is to describe the gravitational field by using the Teukolsky 
formalism |34| implemented in the time domain (using a numerical code introduced in 35]) and by modelling the 
small body by smearing the singularities in the source term by the use of narrow Gaussian distributions. This work 
has also the problem that the size of the small body is too big (in comparison with the size of the MBH) . This shows 
the main underlying difficulty in the numerical simulation of EMRBs, namely, the problem involves a vast range of 
physical scales (spatial and temporal) that expand over several orders of magnitude. Specifically, one needs to handle 
not only large wavelength scales comparable to the massive black hole, but also to resolve the scales in the vicinity of 
the small object where radiation reaction effects play a crucial role. 

An obvious conclusion we can extract from these facts is that, in order to carry out successful numerical simulations 
of EMRBs, we need a high degree of adaptivity. Our proposal is the use of the Finite Element Method (FEM) as 
a natural choice to achieve this high level of adaptivity. Finite Element Methods have not been used occasionally 
in General Relativity (see ,36. 37]) To demonstrate that these numerical techniques have great potential of leading 
to successful simulations of EMRBs, we present results from simulations of a toy problem consisting of a point-like 
source orbiting a black hole in scalar gravitation in 2 -I- 1 dimensions. Our aim is to test FEM techniques in a simple 
representative problem that possesses the main ingredients and challenges of the astrophysical EMRBs. That is, 
at this stage we are basically conducting a feasibility study. In our toy model, the spacetime metric is fixed and 
aims at describing the gravitational field of a non-rotating black hole. For computational efficiency, we have made 
a reduction from three spatial dimensions to two. In this reduction, which we explain in detail later, the metric we 
work with is not longer a solution of the Einstein vacuum equations but it keeps the important property that its 
geodesies coincide with the equatorial geodesies of the Schwarzschild spacetime. This metric is not dynamical but 
fixed. The dynamical gravitational field is described by a scalar field on this spacetime, satisfying a wave-like equation 
with a source term that describes the presence of the small object. An important ingredient of our model is the use 
of a particle description for the small object. The equations of motion of this particle are the geodesies of the fixed 
spacetime metric modified by the presence of (spatial) components of the gradient of the gravitational scalar field. In 
this way we have that the particle orbits the black hole subject to radiation reaction: indeed, the particle generates 
the scalar gravitational field which affects its own motion. 

In this work we compare numerical simulations that use the simple classical FEM with simulations that use an 
adaptive-mesh FEM. The essence of the adaptive mesh technique is to produce real-time local mesh coarsening or 
refinement to achieve the desired level of smoothness in the solution. To that end, a good posteriori error estimator 
to predict the regions in the computational domain where rapidly changes take place is extremely important. There 
are several ways to approach this problem. Theoretical research on this subject [3^ has found that the Hessian 
matrix of the numerical solution can accurately predict where the steepest gradients of the solution would take place. 
We demonstrate that this technique when applied to our toy model easily captures the dynamics of the field in the 
vicinity of the particle since around the particle location either the field or the source term will change much more 
than anywhere else. We then refine the local area surrounding the particle and resolve it more accurately so that we 
are able to achieve our final target. 

The plan of this paper is the following: In Section ^ we introduce the particular theoretical model we want to 
study, namely: the description of the gravitational field and of the point-like source, the computational setup, and an 
energy balance test that can be used to test the numerical computations. In Section lTTll we describe the computational 
techniques that we use for the time-domain simulations of our theoretical model: the FEM, the discretization of our 
equations using Finite Elements, the numerical method for solving the equations of motion of the particle, and finally, 
the Adaptive Finite Element Method ( AFEM) . In Section IIVI we present and discuss the numerical results of the 
simulations. Here, we distinguish between the simulations that use the classical FEM and those that use the AFEM. 
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We devote the last Section to discuss the main resuhs of this work and the future perspectives that it opens. 

II. THEORETICAL ASPECTS OF THE TOY MODEL 

In the following, we describe in detail the main aspects of our toy model, namely: (i) the gravitational theory we 
use, (ii) the derivation of the equations to be solved, (iii) the particle's description, and (iv) the computational setup. 
We use physical units in which c = G = 1 . We use lowercase Greek letters for spacetime indices and lowercase Latin 
letters for spatial indices. In subsections III Al and III Bl indices run as: /i,j^,... = — 3;i,j,... = l — 3. 

A. Scalar Gravity 

Scalar gravitation is a theory of gravity that although it cannot we applied to our physical world (we know it 
cannot fit all the available experimental data), is a good laboratory for numerical relativity due to its simplicity. Of 
particular interest for our purposes is the fact that (scalar) gravitational waves exist in scalar gravity. In this theory, 
the gravitational field is described by a scalar field $(a;'^) on a spacetime geometry described by a non-dynamical 
metric tensor g^^^, which we just prescribe (the version presented in problem 7.1 of the textbook by Misner, Thorne and 
Wheeler ^33] and in pol l4lL \42^ considers only the case in which the non-dynamical spacetime is the flat spacetime) . 
The scalar field does not affect to the spacetime structure defined by g^^ . We are interested in studying the evolution 
of a particle-like object orbiting a black hole in this theory. To that end, we consider a background spacetime metric 
g^y describing the geometry of a non-rotating black hole, and a particle of mass m that follows the worldline z^{t), 
where r denotes the particle proper time. The action of this particle-field system is 

5 = 1 ^Ld\ ^-J^ (^g'^"^V^$V.$ - pe* j d'x , (1) 

where the comoving density p, describing the particle, is given by 

/77? TD 
-^S^ixf^ - z'^{T)]dT = -^^5'[x^ - z^{t)] , (2) 

where i is a time coordinate and u* is the corresponding component of the velocity of the particle: 

Varying the action with respect to the scalar field $ , we obtain the gravitational field equation: 

g^''V^V,$ = 4^e*p. (4) 

Varying now the action with respect to z^ , we obtain the particle's equations of motion 

u'^V^u'' + (g^'' + u^'a'')V^$ = . (5) 

Equations Q and lO) form a coupled system of partial and ordinary differential equations respectively. Equation 
is a hyperbolic and nonlinear equation that describes the dynamical gravitational field whereas lO is a system 
of equations of motion for the particle. In the absence of $ the particle follows the geodesies of the background 
spacetime g^^, but in the present of $ it will not longer follow the geodesic of the background. The way in which $ 
influences the motion of the particle is through its gradient projected orthogonally to the four- velocity of the particle. 
At the same time, the particle influences the gravitational field through the source term p. The coupling between 
the particle and the field is where the radiation reaction is encoded: The particle induces a nonzero $ , and this field 
induces a deviation in the motion of the particle from the background geodesic motion, and at the same this motion 
affects the evolution of the field, and so on. Therefore, the situation is the same as in the general relativistic case 
where the particle is treated as a perturbation on the background spacetime of the MBH, and where the radiation 
reaction mechanism works in the same way. 

This system has a well-defined total energy-momentum tensor , which is given by 
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It has two differentiated components, one associated with the scalar gravitational field, ^*^T^,y, and the other with 
the particle, ^''''T^^: T^^, ~ + ^'''^T^^. Their expressions are respectively 



(7) 
(8) 
(9) 



One can show that the energy-momentum conservation equation 

V^T'^'' = 0, 

follows provided matter conservation holds, that is, V^(pu^) — 0. 

B. 3+1 decomposition of the equations 

In order to numerically solve the equations (|4I5|I it is convenient to rewrite them in a 3+1 language, which allows 
for an initial-value problem formulation. To that end, we follow the usual 3+1 decomposition used in numerical 
relativity. That is, we write our spacetime line-element as follows 



ds^ = -a^df + h^j{dx^ + (3'dt){dx' + 13' dt) 



(10) 



where /ly is the spatial metric of the {t = const.} hypersurfaces, a denotes the lapse function and the shift vector. 
The normal to the hypersurfaces is: 



V = (-a,0), = 1(1,-/3'). 

a 



(11) 



The covariant and contravariant components of the spacetime metric are 

/ 

V 

























hi] I 









(12) 



where h^' is the inverse of hij. Then, h'^y — g^i, + n^ni, is the projector orthogonal to the hypersurfaces {t = const.}, 
which also contains the induced metric that these hypersurfaces inherit: Its spatial components coincide with hij. 

Equation (0J for $ is second-order in space and in time. We are going to split it into two equations that are 
first-order in time. To that end, we introduce the following definition: 



n = = -{dt- P'd,) $ , 

a 

Then, equation Q can be split into the following two equations: 

{dt - I3'd,) $ = an , 

{dt - P'd^) n = a (if n + a'A* + A$ - 47re*p) . 
In the second equation, K denotes the trace of the extrinsic curvature K^^y of the hypersurfaces {t = const.} 



(13) 
(14) 
(15) 

(16) 



The symbols Di and A denote the covariant derivative and Laplacian associated with the induced metric hij. Then 
we can write 



(17) 
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Moreover, a' denotes the spatial components of h^i^n'^Van" . After some calculations we find the following expression 

for a* 

a' = h'^Dj Ina = h'^dj Ina . (18) 

Now, let us perform the 3+1 decomposition of the equations for the trajectory of the particle. Our coordinate 
system {t,x') is, in principle, not adapted to the particle trajectory in the sense that t is not the proper time and 
are not comoving coordinates. Then, in this coordinate system the trajectory of the particle is given by {t, = z^{t)), 
and the unit tangent velocity vector is: 

u^'^|, = u*dt + u'di . (19) 

But also 

U^d^=dr => = = (20) 

dr dr 

Then, on the trajectory of the particle the following holds 

i dz'{t) (it .J , ^ d7j 
= — r^— = v^u* , where = — . (21) 

dt dr dt 

Due to the fact that is a unit timelike vector field {g^^u^u'^ = — 1) we do not need to solve for all the components. 
One possibility is to solve the particle equations for the quantities {z^{t),v^{t)). Then, the equations that we get are: 

j/ = fi+fi, (23) 

where 

/i = (^'r*, - T%) u'^u'^ , (24) 

[U ) 

and the velocity component u* can be written in terms of t;' and the spacetime metric as follows 

(n*)2 = - (g,, + 2g,,^;' + gijv'v^y' . (26) 

The term /| gives the contribution to the geodesic motion in the background spacetime g^^, whereas the term 
describe the deviation from the geodesic motion due to the action of the scalar gravitational field 



C. Pseudo-Schwarzschild Black Hole Background 



We want the background metric g^^ to represent the geometry of a non-rotating black hole. In 3+1 we can use the 
Schwarzschild metric. However, to reduce the computational cost of our simulations, we reduce the space of our toy 
model to two spatial dimensions. There are several ways in which this reduction can be performed. For instance, one 
can start directly with the 2 + 1 metric obtained from the Schwarzschild metric on the hypersurface z = Q. Another 
possibility would be to start from the equations in the 3 + 1 setup, expanding all the different terms and then, to 
neglect all the derivatives with respect to the coordinate z and the components of the different objects in that direction, 
and finally to set z = 0. In performing these dimensional reductions the metrics we obtain are no longer solutions 
of Einstein's equations in the dimensionally-reduced spacetime, but this is not an issue in our toy model since the 
spacetime metric is not a dynamical object. The important thing is that the metric we obtain from these reductions 
keeps most of the important properties of the Schwarzschild metric, in particular, the fact that from the two different 
reductions mentioned above the equations for the geodesies are the same, and they coincide with the geodesies of 
Schwarzschild in the plane z = Q. The equations for when we introduce the expression of the dimensionally-reduced 
metrics, would be in general different. In this work we use the first possibility for the dimensional reduction. As a 
consequence of this reduction the indices we use in this subsection run as follows: jjL ,v , . . . = — 2; i ,j , . . . = 1,2. 
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In what follows we specify the form of the 2+1 background. We start from the Schwarzschild metric in Cartesian 
Kerr-Schild coordinates but reducing a space dimension. That is, our background metric is described by the following 
line element: 

ds^ = rjf.^dx'^dx'' + 2Hif,Udx>'dx'' , (27) 

where 77^^ = diag(— 1,1,1) is the 2+1 Minkowski metric, £^ is a future-directed light-like vector field (both with 
respect to the Minkowski and Schwarzschild metric), and _ff is a scalar. In Cartesian Kerr-Schild coordinates these 
three objects are given by 

X- M 

r]^^dx''dx'' = -dt'^ + S.jdx'dx^ , lf,dx^' = -dt + -j-dx' , H = — , (28) 



where M is the black hole mass and r = ^jK^Jx^ . 

We choose the {t = const.} foliation of the 2+1 spacetimc. Then, the values of the relevant 2+1 quantities that 
result from this choice are: 



a 



2 



1 1 



1 + 1 + 2M ' 



(29) 



= 5,, + 2Hklj , h'' = 6^' - TT^^'^' ' (31) 

i + zn 



Vh = ^Jdet{h,j) = VTT2H = y/ 1 + ^ , (32) 

2M / 2M\"^/^ , , 

K = ^ 1 + , (33) 

1 + 2F" (1+2M)2 ^3 • l-^^i 

As it happens in any Kerr-Schild metric, the determinant of the spacetime metric is equal to —1 , and hence we have 
the following relation between the lapse and the determinant of the spatial metric: \/h = 1/a, which agrees with 
equations (|29I32|I . 

The Dirac delta function that appears in the source p is regularized by using a Gaussian distribution: 

5'[x^^z\t)]^ ^ e-i^, (35) 
(72^ cr) 

where 



^Y.^x^ ^ z\t)f . (36) 



1=1 



The source function p takes the form 



lit, 1 1 L ^ 

p = 2 s ^ — 2 — ^ ^ ■ (3 ' ) 

As we have already mentioned, m is the total rest mass of the particle. We can recover this quantity, which is a 
constant of motion, from the following expression: 



IP 



Vhd^x. (38) 
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where 7 = y/T+lWuiUj = au* . Note that by making the approximation H35|l , the exact character of this relation does 
not change due to the fact that the integral of the Dirac delta and a given Gaussian over the whole two-dimensional 
domain yields the same result. On the other hand, the matter conservation relation V ^{pu^^) — ensures that m, as 
defined by equation (|38l) . is a conserved quantity. 

The next thing wc need is the explicit form of the equations of motions for the particle, which means that we need 
to compute the right-hand side of equation H23|l . or equivalently, the terms (|24I25() . for our spacctimc metric H27|l . 
The expressions that we obtain are: 



M 



2M 
r 



2M 

r 

AM 



+ 6jkV^v'' - 2 + 



M 
r 



2M \ (xjv^) 



(39) 



, ^ ,, 2M I XiV^ 
1 - SjkV^v'^ I 1 



2M 



1 + 2M r r 



2M 



2M 
1 + I W 

r 



n 



1 + ^ 



(40) 



where we have used that 



1 - V^Vi 



2M / x,v' 
r \ r 



(41) 



It is important to remark that the expression (|39|l has the same form as in the 3-1-1 case, confirming the fact that the 
geodesies of our 2-|-l spacetime coincide with the equatorial geodesies (z = 0) of the Schwarzschild spacetime. 



D. Toy Model Setup 

To summarize, we have to solve a set of six equations. Two of them are PDEs, equations (|14I15|I . and contain a 
non linear term, the one that introduces the coupling between the gravitational scalar field and the matter sources. 
The other four equations are ODEs, equations (|22I23|I . describing the trajectory followed by the particle of mass m. 
The two sets of equations are coupled. 

To integrate these equations we are going to consider the following type of spacetime domains: [to,i/] x where 
to and tf are the initial and final integration times, and Vl is the spatial domain (see Fig. the domain at every 
{< = const.} slice, consisting of two circular boundaries: an outer boundary dflout at r = rout > 50M and a inner 
boundary dQm at r = r^n < where — 2M is the horizon radius. The aim is to locate rout far enough close 
to the radiation zone so we can impose standard outgoing radiation boundary conditions. With regard to the inner 
boundary, the idea is to excise the black hole singularity from the computational domain, as it is done in many full 
numerical relativity calculations of black hole dynamics [isHil]. without affecting the computation of the field. This 
can be done as long as we have < rh since the characteristics of the PDEs (|14I15() , for r < rh , all point in the 
direction of r = , which entitles us to perform the singularity excision. Moreover, because of this property of the 
characteristics, we do not need to impose any kind of boundary conditions at the inner boundary. On the outer 
boundary (r = rout) we use the two-dimensional Sommerfcld boundary condition (see, e.g. |45l l46| l: 

= 0. (42) 

We can rewrite this condition, using polar coordinates, like {dt + dr){y/r ^)\r=ro^,t = 0- However, in contrast to 
what happens in the three-dimensional case, where F(t — r)/r is an exact solution coinciding with the radiative 
behaviour of the field (at large r), in the two-dimensional case F{t — r)l^ is not a solution of the equations, and 
therefore the boundary condition (|42|l does not capture correctly the radiative behaviour of the model. In this sense, 
this boundary condition could be improved along the lines shown in references |45l l46j | , by considering higher-order 
derivative boundary conditions, or along the works [4^ l48l l49l l50l | where exact radiative boundary conditions are 
explored. 
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The initial data consists of two functions on R^, ((f>o,9t$o) or equivalently {^0,^0), the initial position of the 
particle zl = z*(io) and its initial velocity vl — v'^{to). For the scalar gravitational field we are going to use the 
simplest initial data, which consists in setting the initial value of $ and its time derivative equal to zero [$o(a;,?/) = 
<^{t = to,x,y), ^o{x,y) = dt<S>{t,x,y)\t=tJ: 

$o(a;,y)=0, and $o(x,y)=0. (43) 

In terms of (<I>, 11) this translates to 

0' 

$o(x,?/) = 0, and Uoix.y) ^ -—d,^oix,y) . (44) 

a 

From a physical point of view, this data corresponds to a situation in which the particle comes into existence at 
the initial time t = to and through the source term in equation H15(l induces a non-zero scalar gravitational field. A 
consequence of this initial setup is the triggering of an spurious burst of radiation. 

The particle's initial data, (z^jti^), is chosen in such a way that it coincides with initial data that, in the absence 
of scalar gravitational field, would correspond to circular geodesies of our background spacetime H27|) . Without loss 
of generality, we assume that the particle is located initially at 

4 = (xo,0), (45) 

where here we can prescribe freely the initial position Xq- Then, the initial velocity is 

(46) 

By using our knowledge of the geodesies [we recall that the geodesies of (|27|l coincide with the equatorial geodesies 
of the 3+1 Schwarzschild spacetime], we can have an estimation of the time that the particle takes for completing a 
circular orbit, which we call T. This time can be obtained by using Kepler's third law: 

47r2 

- , (47) 

where R is the radius of the circular orbit. It is important to remark that this expression is exact when we neglect 
the effect of the scalar gravitational field <&, otherwise it just gives an estimation for the orbital period. For an orbital 
radius r = 10 M we have T « 198.7 M . 
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E. Energy-balance Test 

An important feature of the toy model is the existence of a local conservation law © that involves both the energy- 
momentum of the particle and of the dynamical scalar gravitational field. This local conservation law together with 
the symmetries of the spacetime lead to global conservation laws. We can use these global conservation laws as a test 
for the numerical simulations of our toy model. 

Let us derive the conservation law associated with the timelikc Killing vector field that makes our spacetime static, 
^ — dt- Contracting equation with and using the Killing equations, V^^^ + V^/^^ = , we obtain: 

e^V.T^'-^O ^ V, (T^'-^m) - ?^'"'V,eM = ^ V,(r^%) = 0. (48) 

We are now going to integrate this relation over a compact region of the spacetime, V , and we use the Gauss theorem 
to convert volume integrals into surface integrals: 

f VA^f^Tn dV = =^ f T^-'^^dJl.^O, (49) 
Jv Jav 

where dV is the boundary of V, which is a closed hypersurface, and dS^ the volume 1-form associated with dV. In 
this way we have obtained an energy-momentum conservation law that tell us that the flux of the vector T^^S^'^ across 
the boundary of the region V must vanish. 

The region V is chosen as follows (see Figure 121: It consists of two cylinders, Ci and C2, that extend along the 
timelike Killing direction, and two slices. Si and orthogonal to the timelike Killing, such that they cut the cylinders 
forming a closed 2-1-1 spacetime region. Then, the spacetime region V can be described as 

V = {(x^) ^{t,x,y)\tf>t>U,r2>r>ri,2Tr>0>O}, (50) 

where the angle 9 is defined by tan 9 — y/x; ti and tf can be chosen as the initial and final integration times 
respectively; r2 can be taken to coincide with the outer boundary and ri can be chosen so that ri > 2A/, and hence 
the black hole region is excluded and ^ is timelike everywhere in V. The boundary is then given by 

av = Ci U C2 U 5i U 52 , (51) 

where 

Ci = {{x^') I t/ > t > , r = ri , 27r > > 0} , C2 = {(x^) | t/ > t > , r = r2 , 27r > > 0} , (52) 



5i = {(x^) I t = ,r2 > r > ri ,27r > > 0} , ^2 = {{x^") \ t ^ t f > r > ,2tt > 9 > 0} 



(53) 



When we apply the conservation law (|49|) to our domain, we obtain the following relation between boundary 
integrals: 



J C\ 'J C2 S-\_ ■J S2 



(54) 



This relation can be interpreted by saying that the difference in energy between two given instances of time, ti and 
tf , (which corresponds to the integrals on the spacelike slices Si and ^2) is due to the loss of energy through the 
cylinders Ci (gravitational waves being absorbed by the horizon) and C2 (gravitational waves escaping to infinity). 

In order to evaluate these surface integrals we need first to find the normal vector everywhere on the boundary dV , 
that is, on each of the disjoint pieces of (|51|l . The pieces Si and ^2 [see equation (|53|) ] have constant time t, therefore 
the (timelike) normals there are given by 



dt 



^1 + 2M/r 



dt 



n 



S2 



Si 



yjl + 2M/r 



(55) 
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whereas the pieces Ci and C2 [see equation (I52|) ] are cylinders of constant radius r. Here, it is important to remark 
that both of the cylinders are assumed to be located aX r > Vh — 2M, so that they are timelike hypersurfaces. In 
practice, we are going to take ri very close to . Taking this into account we can write the (spacelike) normals as 
follows: 



jdx^ 



dr 



,dx^ 



dr 



n 



Ci 



n^/l - 2M/ri riv^l - 2M/ri 



n 



C2 



r2 



VI - 2M/r2 r2^l-2M/r2 



(56) 





Si 



FIG. 2: Spacetime region V [see equation (jSOJ] where the conservation law II49|I is tested. 



Then, the contributions of the gravitational scalar field energy-momentum to the surface integrals are: 



o 2Af/r /x^ ^ y ^ 



^1 + 2Af/r Vr " r / J 
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2M/r AT, 
l + 2M/r V7''""^ r"^^ 



4Af/r 



v/1 + 2M/r Vr 



C2 



1 ri 



47r 1 - 


f 2M/ri 




1 




- 2A//ri 
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r2 



dt / de 



n 



2M/ri 



+ 2M/ri 



{cosOd^^ + sin 6195^ $) 



2M^ 



-n 



(cos 6*93;$ + sin9dy^) 



47r 1 + 2Af/r2 
1 



ti -JO 



n 



2M/r2 



^1 + 2M/r2 



2M, 
r2 



-n 



+ 



V^l + 2M/r2 



(cos 09a; <I> + sin 9dy^) 
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The contributions of the particle's energy-momentum to the surface integrals are: 

2M 2M 



Si 



dx dy pe 



I ^P^T^^^^di:, = - I I dxdype 



1 + -^{xy) 



1 + -^{X^V') 



(61) 



(62) 



Since the particle density is very small outside a small neighbourhood around the particle location, the contribution 
due to the particle to the integrals on the cylinders Ci and C2 can be neglected: 



JCl JC2 



0, 



(63) 



where in (|57|l and H61|l . <f> and 11 mean ^{ti, x, y) and Il{ti,x, y) respectively. In (|58|l and H62|l . <i> and 11 mean <i>(t/, x, y) 
and I{{tf^x,y) respectively. Both in H59|) and in (|60|l we have used polar coordinates instead of Cartesian ones (a; — 
r cos 9, y = r sin 6*, dxdy = r dr dO). Therefore, in (|59|l . $ and 11 mean <i>(i, ri cos 0, ri snvO) and I\-{t, ri cos 6, ri ainO) 
respectively; and in H6U|I . $ and 11 mean ^{t, r2 cos 6, r2 sin 9) and Il{t, r2 cos 9, r2 sin 9) respectively. In (|61|) and H62|l. 
the objects p, and it* must be substituted by their expressions l|37|) and H41|l respectively. 
Then, the conservation law that we have to test numerically is: 
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dx dy pe 
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r r- 
2M _^ 2M 



(64) 



In order to decide whether the result from the numerical computation is satisfactory, we can normalize the previous 
expression with respect, for instance, the initial energy of the particle, i.e. the last but one line in this conservation 
law. 



III. THE NUMERICAL FRAMEWORK 



The goal of this section is twofold. First, we want to give a brief introduction to the Finite Element method, a 
numerical technique that has rarely been used in numerical relativistic calculations. A detailed basic exposition of 
the FEM can be found in classical textbooks like 51, 52, 53, 54, 55]. Second, we want to describe the Adaptive Finite 
Element method, where we present a new local mesh refinement technique based on an interpolation error estimate 
introduced recently by Chen, Sun and Xu |38|. This is the technique that we investigate in this paper as a possible 
tool for achieving the adaptivity that the simulations of EMRBs require. 



A. Introduction to the Finite Element Method 



The FEM is a numerical technique for solving problems described by PDEs or that can be formulated as functional 
minimization problems. In what follows, we briefly introduce the basic ideas and procedures of the FEM by using a 
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simple example involving the following wave equation: 



an 



2R, 



^{t^O,x,y)^x^ {dt^)it = 0,x,y) = 0. 



(65) 
(66) 
(67) 



The domain il is a disk of radius R. We prescribe the von Neumann boundary condition H66(l where are the 
components of the normal to the boundary dQ. The initial data is described by H67|) . This problem has the simple 
analytic solution: 



^it,x,y) = t^ + x^ + y^ 



(68) 



We start from the discretization of the computational domain Q into an assembly of disjoint element domains {^2^}, 
that is 



for /? 7 . 



(69) 



In two dimensions the element domains are typically triangles and quadrilaterals. In this work we use only triangles. 
In practice, mesh generation is carried out by using the software included in the general-purpose software package 
FEPG M - A typical mesh (domain discretization) for our example is given in Figure 13 Every element is equipped 
with a finite-dimensional functional space J- a, so that we approximate our physical solution locally, at every element, 
as a linear combination of functions of J^a [usually, a special treatment of the boundary elements is required]. It is 
very common that the functional spaces J-'a are formed by piecewise polynomials. In this paper we consider linear 
elements, where the element functions are first-order polynomials, i.e. a + bx + cy . This choice of element functions 
implies second-order convergence to the solution in the norm. 




FIG. 3: Mesh corresponding to the domain Q in Equation H65|l . 



A very important ingredient of the FEM is that it works with an integral form of the equation we want to solve, 
what is called the weak form of the problem. To obtain it, we multiply the equation (|65|l by a test function 0, integrate 
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over the domain fl, and use the Gauss or divergence theorem which introduces the boundary conditions H66|l . The 
result can be written in the form: 

/:[0,V]EE(0,a2^.) + (V(/., ViA)-2 / r^ds + 2{(j>,l) = 0, (70) 

JdQ 

where (, ) denotes the inner product {u,v) = J^^uvdxdy , V denotes the gradient operator, and s is a coordinate on 
dfl . It is important to note that the third term is the resuh of introducing into the weak formulation the boundary 
conditions H6()l) . 

Using the element functional spaces we can expand our solution in terms of nodal basis functions, nA{x,y) {A — 
1, . . . ,N , being N the number of nodes), which are associated with the nodes or grid points of the mesh. Nodal 
functions use to take the unity value at the node to which they are associated and zero at all the other nodes 
{nAixs) = Sab for any node xb). Since we are going to produce a FEM discretization only in space, the expansion 
of our solution in terms of the nodal functions can be written in the form: 

^PheH^m, Mt.x,y)^Y.'^B{t)nB{x,y), (71) 

B 

where H^{il) is the Sobolev space of functions on fl that are, together with their first and second generalized spatial 
derivatives, square integrable, that is, they belong to L^{il) . We are assuming that iph belongs to H^{rt) for any 
time t € [0 , T] . The subscript h denotes a scale associated with the domain discretization, for instance it may be 
proportional to the square root of the average area of the elements that compose the mesh. In our example h refers 
to the maximum mesh diameter in the whole domain. An obvious property of h, is that the size of the elements goes 
to zero as h goes to zero. 

In a time-dependent problem like the one we are considering, the unknowns are going to be the functions tpA(t) . 
The equations for these functions are obtained from the spatial discretization. In a Galerkin-type formulation of the 
FEM, the discretized equations come from the imposition of the vanishing of the residuals: 

£A=C[nA,^ph]^0, (72) 

which consists in taking (j) — ua in (|70|l . Introducing also the expansion H71|) yields the following system of equations 
for the functions -iJja {t) . 

Mab^^b + ^AB^B = Fa , (73) 

B B 

where the matrix Mab = ("■tIj'^b); the so-called mass matrix, and the matrix Kab = (Vn^, Vn^), the so-called 
stiffness matrix, are symmetric and positive-definite matrices. The vector Fa — —2(nA, 1) -1-2 J^j^ ruAds is sometimes 
called the force vector. Equation H73|l is the outcome of the FEM spatial discretization, which is sometimes called the 
semi-discrete form because it consist of a linear system of second-order ordinary differential equations in time. They 
are usually solved by using Finite-Differences methods. One of the most popular methods for second-order in time 
equations is the Newmark method, which is second-order accurate in time (see, e.g. for details). 

Before we discuss the numerical implementation it is worth mentioning two important features of the FEM: (i) The 
piecewise approximations (piecewise linear in our example) of physical fields on finite elements provide good precision 
even with simple approximating functions. Increasing the number of elements we can achieve the desired precision, 
(ii) The local character of the approximation leads to sparse systems of equations once the problem is discretized. 
This helps considerably to solve problems with a very large number of nodal unknowns. 

The numerical implementation of the equations of this paper has been carried out by using the general-purpose 
software package FEPG 56], which can automatically generate finite element Fortran source code based on component 
programming. It can handle many types of problems, including time-dependent non-linear ones, like the one we are 
interested in. As a test, we have implemented the example described in this section in FEPG (see a snapshot of the 
evolution in FigureQJ. We have also studied the convergence properties of the solution. Here, it is important to point 
out that in an unstructured mesh, like the one we are using in this example, to perform a convergence test is not 
as straightforward as it is for structured meshes. We have to define properly the scale h such that by changing it 
we obtain the correct convergence. Starting from the initial mesh (see Figure ISJ, we call the solution we obtain iph- 
Then, we globally refine this initial mesh by transforming every initial triangular element into four smaller triangular 
elements by connecting the three mid points of each edge. Solving our example equation on this mesh leads to a more 
accurate solution that we call ?Aft,/2, and whose associated scale is h/2. By repeating this refinement process we get 
finer meshes, and by solving our equation on them we obtain more accurate solutions, '4'h/2>'j with associated scale 
. We have checked that the solution we obtain converges quadratically in the scale h to the exact solution Ht)8|) 
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by studying the norms ||V'/i/2'= ~ V'oxactlU^ (see the left of Figure IS)). We have also checked that it convergences 
quadratically in the usual way, without making use of the exact solution, just by comparing the norms — iph/2\\L^ 
and \\iph/2 ~ '^h/iWh'^ (see the right of Figure EJ- 





FIG. 4: Snapshots of the evolution of the solution of the example described by equations Ifci51 - I67j . 



Uh/!-U|l/l|U|,/,-U|| 



Uh/4-u||/||u^/5-u|| 



u^-u^,, / U, ,,,-u. 
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FIG. 5: Convergence properties of the solution of the problem H()5|l - lj67ll . On the left we show the convergence to the exact 
solution 1681 . On the right we show the convergence. 



B. Finite Element Discretizations of the Toy Model 

In this section, we are going to develop a FEM formulation for our toy model in the spirit of the ideas presented above. 
The ingredients of this problem are the following: (i) The computational domain was described in subscction lll DI and 
shown in Figure ^ (ii) The equations that describe our model are the PDEs given by equations (|14I15|I . (iii) The 
boundary conditions and initial data are described in subsection III Dl 

We start from the discretization of the computational domain (sec Fig. We use linear triangular elements. 
The aspect of the resulting triangularization is shown in Fig. El It is worth pointing out that the FEM is specially 
well suited for complex domains. In our case, this allows us to use circular boundaries, which adapt better to the 
characteristics of our problem. This is specially important in the case of the inner boundary, which corresponds to 
the fact that we have excised the black hole singularity from the computational domain. 

The next step in our development is the finite element discretization of the equations H14I15|I . We need to construct 
a weak formulation of these equations. To that end, it is very convenient to rewrite the equations in the following 
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FIG. 6: Triangularization of the spatial domain given in Figure Q 



equivalent form: 

a-^dt<S> ^ a~^P'd^^ +n, (74) 
a-^dtU = a-^P'diU + d^{Vhh'^dj<S>) + a-^a'd^'S? + a-^Kn -^Ana-^e^p, (75) 

where we have used the fact that a\/h = 1 . The main reason for casting the equations in this form is that it brings the 
second-derivative terms into the form of the divergence of a spatial vector, without any additional factors. As we did 
with the example of the previous subsection, we are going to discrctize the system of PDEs (|74I75|) for the unknowns 
$ and n by using a FEM discretization (a Galerkin-type formulation) for the spatial dimensions and by using Finite 
Differences methods in time. In this sense, it is important to take into account that our system of equations is of first 
order with respect to time derivatives and of second order with respect to spatial derivatives. In order to discretize 
these equations, particular attention has to be paid to the convection terms in order to keep them under control. We 
deal with this issue by including additional artificial viscosity into the finite element equations that we obtain. 

For the FEM spatial discretization the finite element space that we use is Sh C H'^{fl), which consists of piecewise 
triangular linear interpolation functions. Then, taking all this into account and operating in a similar way as we did 
in the example of subsection IIII Al the discretized problem that we obtain can be introduced in the following way: 
We need to find ($"+i ,n"+i) e Sh d H^{n) , such that the equations 

(q,-i$"+i, = («-!$", 0) + (At 0) - S{h){At —9,$", At —3,4)) + (Am", 0) , (76) 

a a a 



(a-2n"+i,^) = (a-^jjn^^) + (Ai ^a,n",(^) - 5{h){At^diJr,At^d,^) - At{Vhh'^dj^"+\d,(p) 

+ At (— (^) + At (a-^ifn" 93) - At (47ra~ie*"^' p, if) 
a 

.($»+i_$"+ At)v3ds-Ai/ — -a,$"+Vds, (77) 
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hold for any (0, ip) G Sh and for n = 0, 1, . . . , A'' , NlSi = T , where T = t f — to is the total computational time. At 
is the time-step size we use to evolve $ and 11 . We use the notation (•, •) to denote the inner product in Sh, which is 
defined as {u,v) — J^-^uvdfl (V u, u € Sh)- Finally, to deal with the convection terms, arising from the derivatives 
of the fields along the shift vector, we use a streamline diffusion scheme that includes artificial viscosity and which 
is specially adapted to the FEM. The factor S{h) is the penalty ratio of artificial viscosity, which of course depends 
on the mesh size h. It has to be tuned properly in order to obtain a stable computation and at the same time an 
accurate solution. 

The last terms of ()77|l consist of integrals on the boundaries dQi„ and dQout- The way they come into play is the 
following: In the case of the outer boundary dQ,out, we obtain the last but one integral in H77|l after integration by parts 
of the term with second-order spatial derivatives in (|75|) and imposition of the outgoing boundary condition H42|) . In 
the case of the inner boundary dVLinm we obtain the last integral in (|77|l from the same integration by parts. However, 
the resulting integral does not correspond to a proper boundary condition. As we have already mentioned before, we 
do not need any boundary condition at the inner boundary due the particular structure of the characteristics there 
(remember that the inner boundary is inside the horizon), which all point inwards. Therefore, the resulting integral 
is just an integration of a term proportional to 9^$"+^ on the inner boundary. 

Equations H76I77|) are the basis of the computational procedure that we follow to obtain the solution of our problem. 
The basic algorithm consists in computing first ^n+i from H76|) in terms of $„ and n„. The next step is to introduce 
$,1+1 into the right-hand side of 177|) to obtain n„+i in terms of and n„ together. Then, n„+i has to be used in 
order to compute the value of $ in the next time step, and so on. It is not difficult to show that this algorithm is 
second-order accurate both in space and in time. 



C. Computing the Motion of the Particle 



In the procedure we have just described we have omitted the role of the motion of the particle. It is clear that in 
order to evaluate the right-hand side of (|77|l we need to introduce the position of the particle. And for that, we need 
to solve the set of equations (|22I23|I . These ODEs have to be integrated simultaneously with the PDEs, which means 
that every time we evaluate the right-hand of H77|) we need to evolve them a time step At . And then, to use use the 
outcome of the integration of the ODEs (new position and velocity of the particle) in order to evaluate the source of 
equation [equation (|77|) in the discretized system] . 

The type of numerical algorithm we use to solve the ODEs (|22I23|) has to take into account the particular structure 
of the equations H23|) . which are the non-trivial ones. These equations contain two differentiated terms, /| and f^. 
The first term would give us the geodesic motion around a Schwarzschild Black Hole, and therefore the time scale of 
the changes induced by the term /| is the orbital period of the geodesic that the particle would follow by ignoring the 
radiation reaction effects. The second term contains the gradients of the scalar field <& in the neighbourhood of the 
particle position. These terms are the responsible, in our toy model, of the radiation reaction effects, and therefore, 
the time scale of the changes they induce will be in general much smaller than the orbital time scale. 

Taking into account that these ODEs are nonlinear both in and in , we are just going to use an explicit scheme 
with a time step Atg much smaller than the PDF time step At . That is, we split each PDE time step At in many ODE 
time steps A<s , so that we solve equations (|22I23|I in each time step At^ in order to guarantee the accuracy of the 
solution. This scheme is able to approximate the accuracy that an implicit scheme would have provided. The specific 
discretization algorithm we use in our numerical computations to evolve the particle's position, 2:, and velocity, v, 
from a PDE time step t„, with values (2;", v") , to the next time step t„+i, with values (2;"+^, t>"+^) is given by: 



From m = to m ~ M — 1 : 

<\\ - ^" + At. { /g [<+\ + u , ^1^' , n- , V$"] } , (78) 



,,n+l _ n+1 



where Atg denotes the ODE time step, which is related to the PDE time step by the relation At. = At/M. We 
choose the integer M in order to achieve the accuracy we want. The way in which is computed is by using the 
corresponding FEM piecewise polynomial expansion that follows from the expansion in nodal functions of . 
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D. The Adaptive Finite Element Method 



The AFEM is the apphcation of the classical FEM on a series of local adaptive meshes in order to get more accurate 
numerical results with less computational cost. Starting with a given initial coarse mesh, the adaptive mesh on each 
level is generated locally and adaptively in terms of a posteriori error estimate of the finite element solution. Usually 
local refinement takes place in places where the a posteriori errors are much bigger than elsewhere, or in other words, 
where the finite element solution changes steeply. On the other hand, for those places in which the a posteriori errors 
are sufficiently small, the local derefinement will operate in order to eliminate extra grids because in those places 
the solution changes slowly. Thus, during the finite element computation we can adaptively adjust the mesh density 
without losing the numerical accuracy but reducing considerably the computational cost. 

It is obvious that the key part of the AFEM is the a posteriori error estimate. To have a good such estimate 
means that one can precisely find out in which places the mesh should be refined or derefined without introducing 
much numerical pollution. After this is done, the rest of the procedure is standard. For instance, the mesh bisection 
and the finite element approximation. There is presently a number of works in the literature on these issues (see, 
e.g. 0,0,1^). In what follows we introduce our own a posteriori error estimate, which is an interpolation 
error estimate. Then, we elaborate on our local mesh improvement techniques such as refinement, coarsening, and 
the smoothing strategy, which aim to minimize the interpolation error. 

The interpolation error estimate comes from recent work by Chen, Sun, and Xu j33 |. This estimate can be seen 
as the theoretical foundations of our adaptive mesh techniques, that is, our algorithms are aimed at minimizing (or 
at least reducing) the interpolation error by iteratively modifying the grids. We introduce the estimate through the 
following definition: "Let fl he a, bounded domain in R" . Given a function u £ C^(f2), we say that a symmetric 
positive definite matrix H e M"^" is a majoring Hessian of u if 

|^*(V2u)(x)e| < coe*i?(x)C , (^ e M" , X e 1]) , (79) 

for some positive constant Cq" . Here, ^* denotes the transpose of the vector ^ . We then use the majoring Hessian to 
define a new metric 

Hp^ (detny^H , ip>l). (80) 

There are two conditions for a triangulation T/v, where N is the number of simplexes (generic elements), to be a 
nearly optimal mesh in the sense of minimizing the interpolation error in the norm. The first condition consists 
in asking the mesh to capture the high oscillations of the Hessian metric, that is, H should not change very much on 
each element. This condition can be expressed in a more precise way by means of the following statement: "There 
exists two positive constants and ai such that 

aoC'HrC < (x)e < ai^'Hr^ , (^ G M" , x e f!) , (81) 

where Hr denotes the average of H over r e TJv" . The second condition demands the triangularization TJv to be 
quasi-uniform under the new metric induced by Hp. This condition can also be express in a more precise way through 
the following statement: "There exists two positive constants, Po and such that 

, . max^erlrl 

2 — < Pa (VT e Tat) and — — - < /3i , (82) 

mm^errl 

where |f| denotes the volume of r, and dr.j the length of the J-th edge of t under the new metric Hp. The first 
inequality in H82|l means that each r is shape-regular under the metric Hp. The second inequality means that all 
elements r are of comparable size (also under the new metric), which is a global condition. This means that the 
mesh is more dense at the regions where det_ffp(x) is larger. In j3^], it has been proven the important result that a 
triangulation that satisfies both the local condition (|81l) and the global one H82I) yields a good approximation. This 
result can be expressed in a precise way in the form of a theorem: "Let u be a function belonging to C^(f2), TJv a 
triangularization satisfying the conditions 1^81}} and f^jj)) . and uj the linear finite element interpolation of u based on 
the triangulation T/v . Then, the following error estimate holds: 



^iILp(o) < GN- 



VdetCff) , (83) 



for some constant C — C{n,p, cq, ao, ai, (3o, Pi)- This error estimate is optimal in the sense that for a strictly convex 
(or concave) function, the above inequality holds in a reversed direction." 
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The result expressed by this theorem is the basis of the grid adaptation algorithms that we use in this work. 
Roughly speaking, for a given function u, we will adapt our grids in such a way that the conditions given in (|81|) and 
H82|l are satisfied in the best possible way. One important remark we need to make is that the validity of the theorem 
stated above allows for a few exceptions of the condition H82I) for p < oo (sec for details). This is of particular 
importance since in practice it is very difficult to guarantee that the condition l|82l) is satisfied everywhere. 

The next point is to see how the Hessian matrix of the solution can be obtained when the linear finite element 
approximation is used for the discretization of the PDEs of our problem. Since taking piecewise second derivatives of 
piecewise linear functions will give no approximation to the Hessian matrix, special post-processing techniques need 
to be used in order to obtain a reasonable Hessian matrix approximation from linear finite elements. One of the 
most popular techniques is the patch recovery technique proposed by Zienkiewicz and Zhu (ZZ) |6lL l62j . which is 
based on the least squares fitting on local patches. Results from the application of this technique demonstrate that 
it is robust and efficient. The theoretical reason why the ZZ method works is largely understood to be related to 
the superconvergence phenomenon for second-order elliptic boundary-value problems discretized on a finite element 
grid that has certain local symmetry (see the works of Walhbin [13 and Babuska and Strouboulis ) • These classic 
superconvergence results can be used in order to justify the effectiveness of the ZZ method. A significant improvement 
of this type of analysis was introduced recently by Bank and Xu i65i | . In (64] they give superconvergence estimates 
for piecewise linear finite element approximations on quasi-uniform triangular meshes where most pairs of triangles 
that share a common edge form approximate parallelograms. In l64l| they also analyze a post-processing gradient 
recovery scheme, showing that Q^Vu/i, where Qh is the global projection, is a superconvergent approximation to 
Vu. This result leads to a theoretical justification of the ZZ method for such type of grids (see Xu and Zhang 
for details). 

The gradient recovery algorithm that we use in the numerical examples of this paper is based on a new approach due 
to Bank and Xu |65j |. where they use the smoothing iteration of the multigrid method to develop a post-processing 
gradient recovery scheme. This scheme proves to be very efficient for recovering Hessian matrices. All the meth- 
ods mentioned above can be extended to anisotropic grids with some appropriate modifications, but a theoretical 
justification of such extensions is still lacking. Nevertheless, numerical experiments have given satisfactory results. 

Let us now discuss techniques to improve the mesh quality. We define the mesh quality of a given triangulation T 
in terms of the interpolation error: 

Q(r,M,p) = ||ii--ux.r|liP(s:2) , (1 < P < oo) . (84) 

There are three main ways of improving a mesh: (i) Refinement or coarsening through splitting or merging of 
edges [67ll6^l6pl: fu) edge swapping by replacing sets of elements by other such sets while preserving the position of the 
points f nodes) |7Cll|: and (iii) mesh smoothing, which moves the vertexes of the mesh while keeping the connectivity [tH 
172. .73^ .74] . We derive those techniques by minimizing the interpolation error in the norm. 

For the first method, this can be done by equidistributing the edge lengths with respect to the new metric. Thus, we 
compute the edge lengths with respect to the new metric Hp and mark edges whose length is greater than ri d, where 
ri > 1 is a parameter and d is a fixed edge length, the global average edge length. Then, we connect marked edges 
element-wise according to the different situations that can be given. This is illustrated in Figure The coarsening 
operates like an inverse procedure to the refinement process. It marks the edges whose length is smaller than r2 d, 
where r2 < 1 is another parameter. We then shrink this edge to a point and connect to the vertexes of the patch of 
the edge. 




1 marked edge 2 marked edges 3 marked edges 



FIG. 7: Edge based refinement. 

We consider now the case of edge swapping involving four points {xq,} (a = 1 — 4) constituting two adjacent 
triangles and a convex quadrilateral. Let Ti = A123 U A134 and 7^ = A124 U A234 be two triangularizations, where 
Ao,/37 stands for the triangle made up of the points , xa , and x-^ . Then, we choose the triangulation Ti if and 
only if Q{Ti, u,p) < (3(72, u,p) , for some 1 < p < 00 . In |75|, we show this criteria is equivalent to the empty circle 
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criteria when u(x) = |lx|p . Thus it is an appropriate generation of the edge swapping used in the isotropic case to 
the anisotropic case. 

Finally, local smoothing of the mesh adjusts the location of a vertex in its patch fla, which consists of all simplexes 
containing the vertex x^ , without changing the connectivity. The moving of a vertex to its new location is expected 
to improve the quality of the mesh. Several sweeps throughout the whole mesh can be performed to improve the 
overall mesh quality. By minimizing the interpolation error in fla , we move the vertex Xq to the position x* in such 
a way that 



The derivation of this formula can be found in |75j . 

For the application to our numerical computations we use Qh^u^ and Uh in H85|) . Moreover, we need to make 
choices for the different parameters of this adaptive mesh technique. For the order of the p of the norm we use 
p = 1, the norm. According to [75j |. the norm can catch singularities more efficiently than other norms. For the 
multiple r (ri , r2 in the text) of the global average edge length d under the new metric Hp (we use rd as the threshold 
to determine which edge lengths under the new metric are bigger than it, we will then bisect those edges in the local 
refinement process) we take r = 3 initially. Our local refinement procedure is a nested iteration process. There is 
an outer iteration for which, on each step, we reduce the multiple r till rmm (we take = !)■ The main purpose 
of this iteration is to resolve the singularity as precisely as possible and to control the refinement pollution at the 
lowest level. On each outer iteration step, we perform an inner iteration, and there another number r to control this 
iteration. Usually we do not need this number to establish a criterion to stop the inner iteration of local refinement. 
For the sake of getting the optimal mesh, we let this iteration go until all the edge lengths are smaller than rd as 
before. In practice there are occasions when this may take too much time, and in those cases we set up a criterion as 
the ones before in order to control the number of inner iterations. 



In this section we report the results of the simulations of the toy model. We have distinguished two types of 
simulations. Those that were performed by just using the classical FEM, without adding any extra adaptivity. The 
results are discussed in subsection II V Al And those that were performed by using the AFEM described in the previous 
section. The results of these simulations are discussed in subsection IIVBI 



The mass ratio we have considered is /i = 0.01, that is the particle's mass is m = 0.01 M . The inner boundary is 
located at ri„ = IM, inside the horizon rh = 2M . The outer boundary has been located at rout = 50 Af. Regarding 
the initial data, the only parameters that need to be given, in addition to the ones already given, in order to completely 
specify it [see subsection III D| are the width of the Gaussian that we use to regularize the Dirac delta distribution, 
which we take to be tr = 1 M, and the initial position of the particle, which we take to be {xo, Ho) — (10 M, 0) . Finally, 
we give the resolution that we have used for the simulations. To give a measure of the spatial resolution we describe 
the structure of the mesh. It is composed of iV = 9266 triangles, quasi- uniformly distributed, more concentrated 
near the center and coarsening gradually as we approach the outer boundary. Regarding the resolution in the time 
direction, the step size we use for the time evolution is At = 0.01 M. The algorithm we use to study the evolution 
of the gravitational scalar field $ and the motion of the particle is completely described by the explicit schemes of 
equations (|76I77|I and l(7S|) . 

What we have observed is that the orbit of the particle (remember that the initial data corresponds to a circular 
orbit in the absence of the scalar gravitational field) shrinks gradually until the particle reaches the horizon at 
rji = 2M. The trajectory followed by the particle is drawn in Figure [S] Snapshots characterized by the time t and 
time-step number n of the evolution of the scalar gravitational field are shown by means of the contour plots given in 
Figures [91131 The position of the particle is evident in these graphs. 

We have checked the energy-balance law (|64|) for the parameters given above. The result can be found in Figure IHI 
where the horizontal axis denotes the time and the vertical axis indicates the value of the left-hand side of H64I) in 
units of m . As we can see there, the energy law (I64II is satisfied up to a certain level along the evolution. However, 
from t = iOM the error in the energy-balance test grows and stabilizes around a different but bigger value. This 




(85) 



IV. NUMERICAL RESULTS 



A. Numerical simulations of the toy model using the FEM 



20 



growth is due to outer boundary effects: The particle started from a position r = lOM and the outer boundary is 
located at r = 50 M . As we have mentioned in the discussion of the boundary conditions in subsection lll Ul the outer 
boundary condition l|42|l is not optimal, and induces some error in the solution. This could be improved by either 
moving the outer boundary further out or by using an improved outer boundary condition. 




FIG. 8: Trajectory followed by the particle. 




FIG. 9: Scalar gravitational field $ at time t = 10 M (n=1000) [left] and at time t = 20 M (n=2000) [right]. 

The key point in these computations is the ability of our numerical scheme to resolve the source term describing the 
particle. It enters the equations for the scalar gravitational field $ as a very singular source term. This source term 
depends on the position of the particle, but at the same time, the position of the particle depends of the gradients 
of <i> . Therefore, it is very important first to compute properly the effect of the source p on the field <i>, and then to 
compute accurately the gradient of the field itself. We have already mentioned that in the computations shown in this 
subsection the width of the Gaussian has been taken to be cr = 1 M in order to make the source sufficiently smooth 
for the resolution we have used. In this sense, if we try to use a smaller width, for instance cr = 0.1 M , the classical 
FEM will fail to provide the expected accuracy, although it still provides a numerical solution. Here is where one 
realizes the potential of the AFEM as a better choice to carry numerical computation for a model presenting features 
similar to the ones of our model. In the next subsection we describe how we have implemented the AFEM in the case 
of the toy model and show that it provides a reasonable solution for the case in which ct = 0.1 Af . 
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FIG. 10: Scalar gravitational field ^ at time t = 30M (n=3000) [left] and at time t = 40 M (n=4000) [right]. 




FIG. 11: Scalar gravitational field $ at time t = 50 M (n=5000) [left] and at time t = 100 M (n=10000) [right]. 




FIG. 12: Scalar gravitational field $ at time t = 200M (n=20000) [left] and at time t = 300 M (n=30000) [right]. 




FIG. 13: Scalar gravitational field ^ at time t = 400 M (n=40000) [left] and at time t 



= 500 Af (n=50000) [right]. 
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FIG. 14: Energy-Balance Test for the simulation with the FEM (no additional adaptivity). 



B. Numerical simulations of the toy model using the AFEM 



We have just concluded that the central part of the numerical simulation of our toy model, and more in general 
of EMRBs, is the proper resolution of the matter source, p , and of the gravitational field, $ , in the surroundings of 
the particle position, which is a key issue in order to compute accurately the motion of the particle and hence, the 
waves that are emitted as a consequence. To that end, within the framework we have established above, we want to 
able to compute accurately for small values of the particle's width a. In the limit cr — » we recover the distributional 
description of the particle's that appears in the contimmm description. Then, it is obvious that better accuracy in the 
terms discussed above means to increase the resolution around the particle's position, and in this sense the AFEM is 
a natural choice since it provides the resolution required at the different regions of the computational domain, which 
maintains the computational cost at realistic levels. 

According to the theoretical foundations of the adaptive mesh technique presented above, in order to apply mesh 
adaptivity successfully in our numerical computations, first of all, we have to study which quantity in the problem we 
are dealing with can be used in order to determine where and how the mesh should be refined. In technical terms this 
means that we have to be able to determine the places where refinement (or derefinement) has to take place in terms 
of the Hessian matrix of the selected quantity. We have to look for place where the majoring Hessian matrix is much 
bigger than anywhere else, which are the places where the quantity we have chosen changes very fast. In the case of 
our toy model there is no too much choice. We can either take the field $ or the source p. In the case of it is 
interesting in the sense that it is one of the goals of our computations. However, since it is the solution of a wave-type 
equation, it evolves like a wave in the sense that the profile of the solution will have maxima and minima which will 
be captured by the adaptive mesh procedure, and hence most parts of the domain will be eventually refined, and this 
may not be the most efficient way to proceed. The other choice, the particle's source p , is a better choice for the 
simple reason that we are also solving for the particle's position at the same time that we are solving the PDEs for 
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the field. Therefore, the recognition of the regions that need refinement is going to be simpler in our problem. That 
is, the most efficient way of adapting the mesh for our computations is to concentrate on the region where the particle 
is present and follow it from there. 

If we look at the behaviour found in the previous subsections (see Figures I^I1U|) . it happens that the bigger values of 
<i> occur around the particle's position. Therefore, refining there is very convenient at that stage of the computation. 
A possible drawback of using the particle's source as the refinement criterion is that the situation changes at later 
times when the field $ is present all over the computational domain and it may be a need for refinement in other 
places like close to the horizon. However, in the numerical simulations we have performed we find that using the 
source term p for refinement is much better than using the field $ . By calculating the Hessian matrix of the particle's 
source we can exactly know where and how to get local mesh refinement on the particle. 

For the numerical simulations of the toy model with the AFEM we use the same parameters as in subsection II V Al 
excepting for the width of the Gaussian, which we have been able to reduce down to cr = 0.1 M . And this has been 
done using a number of triangular elements in the interval 11000 — 15000 (the number of elements changes with the 
evolution, depending on how the adaptivity is implemented), a number comparable to the one used in the calculations 
without adaptivity. That is, working with the AFEM we have been able to use a width a an order of magnitude 
smaller than with the classical FEM, using a comparable number of elements. The time step is now At = 0.005M, 
half the one used in the previous simulations. The trajectory followed by the particle is drawn in Figure [TBI We show 
the evolution of the scalar gravitational field in Figures ^1201 These are contour plots as the ones we used previously. 
In order to show how we refine the mesh in the evolution we have superposed the mesh to the contour plots. In this 
way we can see how the high-resolution part of the mesh moves with the particle. One can clearly see the difference 
between this case and the previous case: the particle inspirals much faster in this case, and it hardly completes more 
than one orbit. To understand this difference, it is important to remark the fact that the particle sources are different 
since the Gaussian profiles have a different width a and the fact that the sources couple non-linearly with the scalar 
gravitational field. As a consequence, the difference in a leads to a different evolution. We have checked that when 
take a = IM with the AFEM we recover the type of trajectories we get in the case without adaptivity. 



A more detailed view of the structure of the refined region around the particle is shown in Figure 221 where one 
can see the different levels of refinement, and how the resolution increases as one approaches the particle. A global 
perspective of the structure of the mesh, containing the excised region around the singularity of the black-hole, and 
the refined region around the particle can be seen in the three-dimensional graph shown in Figure 1^ 

We have also checked the energy-balance law (|64|l for the simulations with adaptivity. We have plotted the result 
in Figure 1^ where again the horizontal axis denotes time and the vertical axis the value of the left-hand side of H64() 
normalized with the mass m of the particle. The behaviour of the error in the conservation law is very similar to the 
one in the case without adaptivity, where the error in the energy-balance test grows after the boundary affects the 
particle. In this case, the error does not seem to stabilize as clearly as in the previous case, and in average the error 
in the conservation law is bigger. However, one has to take into account that the evolution in this case is much faster 
since the particle has plunged very quickly. Moreover, since the number of elements in both cases is comparable, 
but in the case of the AFEM a considerable amount of them have been used to increase the resolution around the 
particle, this means that the density of elements near the boundaries is less and therefore if errors are propagated 
from there, they can affect more the result than in the case without adaptivity. This can be avoided by using more 




FIG. 15: Trajectory followed by the particle. 
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FIG. 16: Scalar gravitational field $ at time t = 5M (n=1000) [left] and at time t = 20M (n=4000) [right]. 




FIG. 17: Scalar gravitational field $ at time t = 35 M (n=7000) [left] and at time t = 50 M (n^lOOOO) [right]. 




FIG. 18: Scalar gravitational field <1> at time t = 65 M (n=13000) [left] and at time t = 80 M (n=16000) [right]. 

computational power to push the boundaries far enough so that we can neglect their effect during the significant part 
of the evolution. 



V. REMARKS AND CONCLUSIONS 

In this work we have studied the application of Finite Element techniques to the time-domain numerical simulations 
of EMRBs. More specifically, we have shown how Finite Elements can help us in achieving the degree of adaptivity 
that the computation of radiation reaction effects around the small body requires. However, adaptivity is not the 
only advantage in using Finite Elements for the study of EMRBs, some of which are obvious from the study of the toy 
model we have considered in this work. In short, the main additional advantages that Finite Elements provide are: 
(i) Versatility, the FEM can be applied to a wide range of problems: static, quasi-static, transient, highly dynamical, 
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FIG. 19: Scalar gravitational field $ at time t = 95 M (n=19000) [left] and at time t = llOM (n=22000) [right]. 




FIG. 20: Scalar gravitational field $ at time t = 125 M (n=25000) [left] and at time t = UO M (n=30000) [right]. 




FIG. 21: Scalar gravitational field $ at time t = 150 (n=25000) [left] and at time t = 160 Af (n=32000) [right]. 

linear and nonlinear, etc. Moreover, the practical implementation of the FEM can be easily designed in a modular 
way, as it shown by the existence of a number of multi-purpose Finite Element packages. This is a good property 
for the design of complex numerical codes, as the ones one may need for the description of astrophysical EMRBs. 
(ii) Many of the procedures that one uses in the framework of the FEM have a solid theoretical foundation, in the 
sense that behind them one often finds rigorous mathematical analysis (examples of this can be found in the classical 
textbooks |5lll53 . l53Ll53 . l55l l). (iii) The ability of the FEM to manage problems with complex geometries. For instance 
we can easily consider spherical boundaries which is helpful in imposing boundary conditions and also in implementing 
excision techniques for black hole singularities. This is something that may be worth to export to fully numerical 
relativistic calculations of spacetimes containing black holes, (iv) In relation with the previous issue, the imposition 
of boundary conditions can be done in a somewhat natural way in a FEM framework. The idea is that by adapting 
the mesh to the geometry of the problem under study, boundary conditions which are also adapted to that geometry 
will be incorporated in a natural way into the FEM discretization through the weak formulation of the problem. The 
outgoing boundary conditions that we have used in our toy model are a good example, its implementation is simpler 
that it would have been in a similar Finite Differences scheme, (v) Another advantage of the FEM, that we have not 
explored in this paper, is that it is to some extent natural to handle distributions, like for instance the typical Dirac 
delta distributions that appear in physical systems containing particles (e.g., EMRBs). The idea is that we can use the 
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FIG. 22: Details of the refinement used around the particle's position. 




FIG. 23: A three-dimensional visualization of the mesh including the excised region and the refinement around the particle. 

exact properties of the distributions at the level of the weak formulation of the equations, which in the case of Dirac 
delta distributions leads automatically to a regularization of the singularity associated with that kind of distributions. 
We can apply this idea, for instance, to the calculation of the general-relativistic perturbations of a non-rotating black 
hole, produced by the presence of a point-like object, via the Regge- Wheeler- Zerilli formalism. Then, the type of FEM 
discretization we would get is similar to the one obtained by Price and Lousto [t^ [t^, Izll , where they also used an 
integral form of the equations in order to obtain a discretization. In addition, the FEM procedure can be generalized 
to higher-dimensional problems in a straightforward way. These ideas will be the subject of future investigations 

In this paper we have illustrated the capabilities of the FEM and its suitability to simulate EMRBs, by applying 
it to a simplified model, a toy model that retains the basic ingredients and difficulties of general relativistic EMRBs. 
We have shown simulations of this system based on the classical FEM and simulations based on the AFEM. An 
outcome of this work is the realization that in order to increase the resolution around the small body, treated as 
a particle, the introduction of adaptivity in the region of the particle is necessary and in this sense the AFEM is 
a natural tool to use. The primary benefit of the AFEM is that it can provide an efficient, accurate and reliable 
computational analysis of very large continuum problems, for only a relatively small fraction of the cost associated 
with the non-adaptive FEM. The accuracy of a finite element solution is directly dependent on the number of free 
parameters used to mathematically represent the problem, and how effectively those parameters, or mathematical 
degrees of freedom, are distributed over the problem space. Furthermore, the full computational cost associated 
with obtaining a finite element solution is related to both the number and the interconnectivity of the degrees of 
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FIG. 24: Energy-Balance Test for the simulation that uses the AFEM. 

freedom used in the problem discretization. Consequently, the most efficient distribution of degrees of freedom for a 
problem is the one that provides a sufficiently accurate solution for the lowest number of free parameters. Currently, 
the only practical way to achieve this objective is by using adaptive solution strategies which are capable of cleverly 
evolving and improving an efficient distribution of degrees of freedom over the problem domain by establishing solution 
error distributions, and then adjusting or adding degrees of freedom to the discretization in order to correct them. 
By increasing the numbers of degrees of freedom only in the regions of higher error in the solution, it is possible 
to make the most significant improvement in the global accuracy of the finite element solution for the minimum 
additional computational cost. In this paper we employed the Hessian matrix of the solution to describe the solution 
error distribution, which can perfectly guide us to where and how we have to locally refine the mesh without any 
unnecessary pollution. In this sense, it is important to emphasize that the rest of adaptive strategies available do not 
have this advantageous property. Moreover, this comes without paying an extra price since the main advantages of 
the other adaptivity techniques are present in the AFEM that we have used. 

The natural continuation of this work is the transfer of the technology used here to the general rclativistic problem. 
One possible way is to consider the general framework of metric perturbations in the setup of the 3+1 decomposition, 
that is, to use 3D perturbative numerical relativity, trying at the same time to profit from the experience gained in 3D 
full numerical relativity. However, 3D calculations using the AFEM may be at present computationally too demanding 
and therefore, other avenues should be also explored. Among them we can consider ID calculations restricted to the 
case of a non-rotating MBH, by using the well-known Regge- Wheeler and Zerilli-Moncricf formalisms, or just by using 
a harmonic gauge, where the computation of self-forces seems more natural. This would be an interesting benchmark 
to test further the FEM techniques in a general-relativistic context. From here, one can study problems of more 
interest for gravitational wave physics related to LISA (involving spinning MBHs) by considering 2D calculations 
using the curvature based formalism of Teukolsky for linear perturbations around Kerr black holes. One can go 
beyond by considering also perturbations of Kerr black holes using metric perturbations in a harmonic gauge, where 
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the computation of self-forces can be carried out by using techniques already present in the literature. 
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